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ABSTRACT 


A three dimensional transient computational model of heat transfer during gas 
tungsten arc welding is generalized, and then validated by comparison to Rosenthal’s 
solution for moving point sources of heat. The current version of the code allows much 
greater flexibility in the specification of the thermal input from the arc. The resulting 
surface temperature profiles and fusion zone shapes are compared to those measured 
experimentally for several input power levels for autogenous gas tungsten arc welding. 
Arc efficiency is experimentally determined using change of phase of a liquid 
fluorocarbon. The model is shown to be useful for modeling autogenous welding of 
thick plates. Weld seam misalignment and surface flaw detection are shown to be pos- 
sible ahead of the arc with accurate surface temperature detection methods. The po- 
tential of the model for creating a database of fusion and heat affected zone sizes, 


temperature profiles, and cooling rates for various materials, processes, and power levels 


is indicated. 
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I. INTRODUCTION 


The pursuit of an automatic adaptive control system for the gas tungsten arc 
welding (GTAW) process has been ongoing for nearly two decades. Numerous re- 
searchers have assisted in developing the state of the art during this time. However, the 
pursuit continues, as many of the so-called automatic welders are not truly automatic; 
- an experienced, highly-skilled operator is the key to the success of these “automatic” 
welders. 

To design a completely automatic system, the first step is to model the process as 
accurately as possible. To obtain an accurate model, all aspects of the process must be 
fully understood. Once the process is understood, mathematical formulas can be devel- 
oped expressing the relationships between the variables in the process. Perhaps the first 
and the most well known of the many models is that of Rosenthal [Ref. 1]. His model 
of the process is the analytic solution of the heat diffusion equation in the presence of 
point, line, and plane moving sources of heat. More recent work has shown that while 
Rosenthal’s solutions are reasonable approximations of reality far away from the heat 
sources, in the vicinity of the fusion zone these are in considerable error. Effects of 
melting and the resulting variations in the thermophysical and transport properties must 
be accounted for in addition to a heat source of finite size in order to estimate the time 
temperature histories. 

The advent of modern control theory, nearly a decade ago, created the opportunity 
to advance the modeling of GITAW and gas metal arc (GMA) welding processes. 
Moody was perhaps the first to publish models of the GMA and GTA welding processes 
in the state-space. He presented descriptions of the process variables and the resulting 
nine first order non-linear differentia! equations relating the variables. His work appears 
to be the starting point for a fully automatic system. One only has to determine the 
exact nature of the relationships to create the system. The effort to determine these re- 
lationships sums up the next decade or so of research in the area of automatic control 
of welding. Not all of the variables in Moody’s model are measurable, which is not a 
problem as long as the unmeasurable variables can be estimated through the use of 
empirical relationships, or some other means. Once all variables are measured or esti- 
mated, a closed loop control system may be designed using the relationships between the 


variables and their derivatives. [Ref. 2] 


A refinement of Rosenthal’s solution, using all but one of his assumptions was pre- 
sented by Eagar and Tsai in 1983. The assumption not used was the point source of 
heat. Eagar and Tsai used a distributed source of heat with a Gaussian distribution. 
The comparison of their results with experimental data revealed a closer agreement than 
Rosenthal’s analysis. They concluded that the strength of the theory presented was that 
it gave accurate functional relationships between process parameters and materials pa- 
rameters. [Ref. 3] 

At nearly the same time, Goldak er al. were developing a finite element model for 
welding. In their paper they take the concept of a distributed source a few steps further 
by applying it over a volume of the work piece to account for the arc digging and stir- 
ring. They note the limitations of applying a disk source at the surface of the material 
to a deeply penetrating process such as electron beam welding. Thus the hemispherical 
power density distribution model is proposed. Then, noting that few if any actual 
molten pools are hemispherical, they proposed an ellipsoidal power density distribution. 
This idea was further generalized to a double ellipsoidal power density distribution. This 
came about because the asymmetry of the temperature gradients ahead and behind the 
arc were not closely matched by the ellipsoidal model. The double ellipsoidal model al- 
lows for using quadrants of two different ellipsoids, front and rear, to more closely model 
actual results. The power of this model for modeling seam welding of two dissimilar 
metals bv splitting the quadrants into octants down the arc centerline was 
noted. [Ref. 4] 

Kou and Sun studied the fluid flow during stationary welding, to later apply their 
results to moving sources of heat. The driving forces for fluid flow in the weld pool for 
aluminum and steels were listed in order of their relative effect on the flow pattern. 
Their work could provide substantial evidence for use in determining the dimensions of 
a double-ellipsoidal power distribution region for a given set of process parameters once 
their results are applied to moving arcs. [Ref. 5] 

_ Tsai and Eagar related current, arc length, electrode tip angle, and shielding gas 
composition to heat input magnitude and heat source distribution. The arc length was 
shown to be the most effective in shaping the heat distribution and the current was 
shown to dominate the magnitude of the heat flux on the surface. They also determined 
that the heat flux was closely approximated by a Gaussian distribution and gave the 
approximate range of the distribution half-width as arc length was increased. [Ref. 6 ] 


Kou and Wang achieved excellent agreement between predicted and measured weld 


pool shapes using a three-dimensional computer simulation of convection in the weld 





pool. The only drawback to their work was that symmetry was invoked, which removes 
the ability to model asymmetrical problems associated with welding such as weld seam- 
electrode misalignment. [Ref. 7] 

The effects of changes in machine variables were investigated at length by Giedt. 
The wide range of efficiencies of the GTA process was explained as due to the fact that 
conduction heat transfer models do not take into account the convective heat transfer 
in the weld pool. Giedt also listed some computing time requirements of interest. Using 
TACO3D, a three-dimensional finite element conduction heat transfer code, with con- 
stant average thermal properties, one hour of CRAY computer time was used. This was 
increased to nearly ten hours using variable properties and a smaller element size to 
smooth contours. One of Giedt’s conclusions is that the effect of weld pool convection 
cannot be accurately represented in a pure conduction model. He notes that use of a 
fictitiously higher thermal conductivity (2 to 5 times that of the solid) has been used to 
account for the effect of convection, but the applicability of doing so is dependent on 
the direction of flow in the weld pool, which may not be known. [Ref. 8] 

Oreper and Szekely examined the development of the weld pool for a tungsten inert 
gas (TIG) spot weld. They noted that during the initial period of the arc, convection 
played very little part, and became more significant as time went on. This behavior was 
not common to all materials. For titanium, neglecting convection in the weld pool ini- 
tially would cause large errors. Three guidelines were presented: 


* convection is not likely to be important for good conductors and shallow weld 
pools. 


® convection is likely to be important for poor conductors and deep weld pools 
* surface tension effects are likely to be important in affecting the circulation pattern 
and subsequent weld pool shape. 
The conclusions stated were that the conditions for which convection plays an important 
part can be defined and that the relative importance of surface tension, buoyancy and 
electromagnetic forces in determining the weld pool shape and circulation patterns can 
be determined. [Ref. 9] 

Lu and Kou used a technique called non-parametric minimization to determine the 
power and current distributions in a gas tungsten arc. This technique was applied to the 
spht anode method. They noted that the power- and current-density distributions in gas 
tungsten arcs are generally steeper than the Gaussian distribution, making it a poor ap- 


proximation in the case of a larger arc gap. The Gaussian distribution was fairly close 


for the smaller arc gaps. [Ref. 10] 





Zacharia et al. developed a computational model code called WELDER. The code 
has nine special features which allow it to realistically simulate many interesting prob- 
lems including inclined welding and welding in microgravity environments. The code 
takes into account the following physical phenomena: 


e melt surface is deformable--the weld crown, surface mnpples and surface gravity 
wave phenomena can be simulated. 


e the local transient effect of phase transformation is incorporated into the energy 
equation. 


e marked-element technique--accurately simulates transient development of the 
solid-liquid interface. 


e incorporates effect of surface tension due to local curvature conditions of the 
deformable surface of the molten metal. 


® considers the magnetohydrodynamic force terms in the momentum equation. 
e ability for motion of the arc and arbitrary positions (inclined welding). 

e different gravitational forces. 

© arbitrary geometry and finite thickness. 


© does not require use of the Boussinesq approximation. 


The capabilities of the code appear to be extensive. CPU times of about 30 minutes were 
reported, but not correlated to the amount of real time simulated. Since the sample size 
(24mm by 24m by 6mm) and material (Al 6061) properties may have had a significant 
effect on the CPU time required, obtaining timely results for a wide range of processes 
and materials is uncertain. [Ref. 11] 

Zacharia er al. reported on the influence of surface active agents and the temper- 
ature distribution on the weld pool surface. 1t was concluded that penetration and the 
weld pool aspect ratio were dependent on a combination of the two factors. The sulfur 
contents (90 ppm and 240 ppm) had a significant effect on the weld penetration for the 
GTA process for 304 stainless steel, however, there was no noticeable difference for the 
laser beam process which was also investigated. [Ref. 12] 

A three-dimensional numerical model of the Gas Tungsten Arc (GTA) welding 
process was created by Ule [Ref. 13], as a tool for the development of an expert welding 
system. The David Taylor Research Laboratory is developing such a system for welding 
submarine hulls. Before using Lle’s model in any application, extensive testing and 
validation was required. 


Testing the numerical model’s Fortran codes was accomplished in several steps. 


These included: 





Comparing the expressions in the Fortran codes to the governing equations. 


Running the programs. 


Comparing the output to previously published output. 


The validation of the model also consisted of several steps. These were: 
e Determination of the arc efficiency of the GTA process. 


e Comparison of the surface temperature profiles with the profiles from the 
Rosenthal solution with the same heat input. 


e Comparison of the surface temperature profiles with the instrumented plate exper- 
imental results for the same heat input. 


e Comparison of weld pool shapes during start-up from experimental results with 
shapes predicted by the codes. 

The model was shown to be useful in predicting the temperature distribution during 

GTA welding. There is continuing research to verify cooling rates during welding 

start-up and shutdown and to develop a low-cost surface temperature sensing system. 


Completion of these research efforts will facilitate the design and subsequent use of an 


expert welding system for GTA welding which with some effort could be adapted to the 
Gas Metal Arc (GMA) welding process. 





Il. MODIFICATIONS TO THE MODEL CODES 


A. MATERIAL AND GRID PARAMETERS 

The first modification to Ule’s basic welding program [Ref. 13: pp. 84-92], was to 
replace parameters which had been hand calculated and entered as constants, with vari- 
able expressions, thus allowing changes in material and welding parameters to be more 
easily effected. Examples of the parameters that were changed in the way they were re- 
presented in the code are: Fo(1), Fo(2), Fo(3), Bi(1), Bi(2), Bi(3), Bi(4). To make changes 
to these variables required the introduction of eighteen additional variables in the pro- 
gram. Since these variables were involved in calculations only once, the additional cost 
ia computation time was negligible. The changes in the Fo(1) expression demonstrate 
how the code was modified for the variables listed above. The Fo(1) expression in Ule’s 


program WELD was: 
Fo(1) = 0.1636Az 


The revised expression used in the START series of programs and all subsequent revised 
codes used the following expression: 


Fo(1) = ——24L__ } 
AxAr(.001) 
where: 
4 is the diffusivity 
At is the time step 
Ax, is the control volume length in the x-direction in the coarse zone 
Ay. is the control volume length in the y-direction in the coarse zone 


(.001) converts from millimeters to meters 
Similar changes were made to the other Fourier and Biot numbers. 


B. GAUSSIAN DISTRIBUTED HEMISPHERICAL HEAT INPUT 
The most significant change to the codes was the manner in which the heat was in- 


put to the metal. Ule had used a Gaussian distribution over a rectangular parallelepiped, 
though the volume was stated to be hemispherical [Ref. 13: pp. 14,86]. The START2 
code was a version of WELD modified to simulate a hemispherical heat input of 4.5 





millimeter radius, with a Gaussian distribution. The Gaussian distribution was modeled 
using the approach outlined by Goldak er al. [Ref. 4: pp. 301-302] for a finite element 
model. The discretized equation for heat input per unit volume became: 


Hit = Hi, + | Ge 
AK Lk 
Gist 
where 
n is the current time level 
i.j,k are the nodal indices 
Gin is the Gaussian factor, e-ME?+(FP +(e 
Vita is the nodal volume 
Q is the power input in watts 
Bij. is the enthalpy per unit volume at a particular time 


It was noted that the constant in the exponent of the Gaussian factor was 1.7 in Lle’s 
codes [Ref. 13: p. 86] which corresponds to a relative intensity of the arc power of about 
18.3 per cent at the edge of the arc, vice five per cent as suggested by Goldak ez al. [Ref. 
4: p. 301]. Using three as the constant gives a five per cent value at the arc boundary. 
An investigation as to the effect of the exponent on the surface temperature profiles was 
made. The exponent was varied between 1.0 and 3.0 in increments of 0.5 and the re- 
sultant surface temperature profiles plotted and visually compared. It was noted that 
the size and shape of the weld pool (near the 1750 K isotherm) was virtually unchanged. 
The differences between the plots were the temperature gradients and peak temperatures 
within the weld pools. As expected, both were higher for the more concentrated heat 
input (using three as the constant). A collection of the plots is provided in Figure 1. 
The variable GAUSS refers to the constant in the exponent. 

Fourteen more variables were added to code this equation. The values of 
a, b,and ¢ are the x,y, and z radii of the arc power distribution in millimeters, respec- 
tively for a continuous model. Discretizing the model enlarges the radii in the 
x- and y-directions and decreases the radius in the z-direction by one-half the nodal 
spacing, or 0.5 millimeters in this case. The limits of summation for the above equation 
were found by dividing a, 6, and c by the node spacing and adding and subtracting this 


value plus one to the arc center. The one was needed to ensure that no nodes were left 


out since the arc position in the y-direction was discretized by the Fortran function INT. 





Adding the one also allows for values of the arc radius other than those ending in .5 and 
for changes in the x-position of the arc as occur in the misalignment programs. As an 


example, the limits of summation for a 4.5 millimeter radius hemispherical heat input 


were: 
x from XARC—5 to XARC+5 
y from YARC—5 to YARC+5 


from 1 to 5 


C. LACK OF FUSION PROGRAMS 

The lack of fusion program, WELDLF, was first altered to input the heat of the arc 
in a hemispherical fashion identical to the START series. The temperature deviations 
from quasi-steady conditions oscillated excessively near the flaw zone. Diagnostic test- 
ing revealed that the variable LOF, used to control the position of the lack of fusion 


zone, was incorrectly determined. The equation used by Ule [Ref. 13: p. 33] was 
LOF = 65—INT(WVEL x TIME) 


which moves the flaw in one millimeter increments as the distance changes. The error 
arises from the fact that the arc is also moving toward the flaw. To correct this, the flaw 
was moved only when the fine zone was shifted, thus ensuring the flaw maintained a 
fixed position in the metal. The equation used in WELDLI and subsequent revisions 


used a variable already defined. 
LOF = 67—STEP 


The variable STEP is incremented by three each time the fine zone is shifted. The dif- 
ference between the two methods of placing the flaw in the material is illustrated by the 
plots in Figure 2. The YARC—JAJIN curves represent the distance from the arc to the 
node at the front of the flaw; the FARC—JAAX curves the distance from the arc to the 
rear of the flaw. The front of the flaw is the edge nearest the arc when the flaw is in- 
troduced into the material. The distance between the curves on each plot corresponds 
to the distance between the first and last nodes of the flaw. The actual flaw distance is 


one node spacing longer, or 4:mm in this case. The effect of the change in the LOF var- 


iable on surface differential temperature profiles is demonstrated by the plots in 
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Figure 1. Surface temperature profiles varying the constant in the Gaussian 
factor; a) GAUSS = 1.0, b) GAUSS = 1.5, c) GAUSS = 2.0, d) 
GAUSS = 2.5, e) GAUSS = 3.0 
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Figure 2. Arc-to-flaw distances: (a) incorrect LOF variable (b) correct LOF 
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Figure 3. This change removed much of the jaggedness of the plotted differential sur- 
face temperatures near the flaw. 

An additional change to the WELDL series was the thermal conductivity of the 
flaw. No material could be found with an extremely low conductivity for which thermal 
conductivity and heat capacity versus temperature relationships over the range of tem- 
peratures encountered during the welding process were available. The conductivity was 
therefore set arbitrarily to a constant 5.0 W/m- K and the heat capacity to the same re- 
lationship as the parent material. This was done in order to examine the effect of a low 
thermal conductivity inclusion on the temperature variation within the material. 

A high conductivity flaw material was also utilized, in a different version of the 
WELDL code series. The material chosen was tungsten. Piecewise linear relationships 
were fit to data from a table in Incropera and DeWitt [Ref. 13] for the thermal 
conductivity and specific heat of tungsten. While the thermal properties of the flaws 
encountered in applications may vary, the intent in this study was to determine the re- 
sponse of the temperature patterns to a flaw of thermal conductivity higher than the 
parent material. Plotted results comparing the low and high conductivity flaws are in- 
cluded in Figure 4. 


D. SURFACE HEAT INPUT 

The program STARSU was utilized to obtain the results of a surface heat input. 
In this case, the heat was input only to the surface nodes in a Gaussian distribution over 
circle of radius 4.5mm. The depth radius was set to zero and the code altered to prevent 
division by zero. 

The results of the surface heat input were unrealistic. The maximum temperatures 
experienced in the weld pool were near 4000 K. One interesting point was that the size 
and shape of the weld pool boundary was relatively unchanged from that obtained for 
the hemispherical heat input. 


E. DOUBLE ELLIPTICAL HEAT INPUT 

A version of the START program, STARDE, was used to observe the effect of a 
double elliptical heat distribution. The program was altered to include the discretized 
form of the equation presented by Goldak er al. [Ref. 4: p. 302]. The double ellipsoid 
used had a 9.0mm minor axis (x-directions), semi-major axes of 4.5mm to the front 
(+ y-direction) and 8.5mm to the rear (—y-direction) and a 1.5mm depth (z-direction). 


The heat was input using a Gaussian distribution with 70 per cent of the net input in the 


rear portion of the ellipsoid and 30 per cent in the front portion. The size and shape of 
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the distribution as well as the front and rear fractions were arbitrary, but not without 
knowledge of the values chosen in other programs for the laser welding process. 

The plotted results were compared to the rectangular parallelepiped, hemispherical, 
and surface heat inputs. It was noted that the size and shape of the fusion zone on the 
surface was not changed much between the surface and hemispherical heat inputs. The 
double elliptical heat input resulted in more oblong weld pools than the other methods. 
The ratio of the semi-major axis ahead of the arc to that behind the arc may have been 
too high as the temperature gradients were not as steep in front of the arc as would be 
expected. The cross sections of the fusion zones were more noticeably different. The 
hemispherical heat input resulted in profiles which appeared deeper than the half-width 
of the fusion zone at the surface. The double elliptical heat input resulted in fusion 
zones with depths less than the half-width at the surface, which was found to be close 
to the shapes observed in later experiments. 


F. ROSENTHAL VERIFICATION CODE 
The program STARTR was created from the START series to simulate the 
Rosenthal solution [Ref. I: pp. 849-869]. This was accomplished in the following man- 
ner: 
¢ reducing the heat input to prevent melting 
® using constant thermal properties (conductivity and specific heat) 
e eliminating convection and radiation from the surfaces 


* adding the heat to a single node on the surface 


The Rosenthal solution was coded in the program ROSEN and the output compared to 
the STARTR output using plots from a revised plotting program. These plots appear 
in Figure 5. 

Rosenthal surface temperature profiles which would have included the infinite point 
(the center of the arc) were drawn using a point-to-point curve drawing routine with the 
infinite point replaced by a large finite value, while all other plots used a natural spline 
curve drawing routine. 

The model adequately predicted the temperatures at the four locations chosen. The 
errors Were maximum at the center of the arc where the Rosenthal solution is known to 
be lacking in its representation of actual conditions. Other than the center of the arc, 


the maximum difference between the model results and the analytic solution was about 


one degree which translates to less than one per cent. 
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G. COOLING RATE CODE 

The program created by Ule [Ref. 12: pp. 127-134] to determine the cooling rates 
during the welding process was modified by changing the heat input to a Gaussian dis- 
tributed hemisphere as was done for the previous codes. Giedt [Ref. 8: pp. 35, 39-40] 
had pointed out that modeling the process with a pure conduction model was inaccurate, 
but that some reasonable results had been obtained in a conduction model by adjusting 
the thermal conductivity of the liquid metal. It was also pointed out that the weld pool 
circulation pattern must be known to properly account for convection when using a 
conduction model. Taking these points into account, a new version of the code, 
WELDC1, was written. The heat was input on the surface in a 4.5 mm radius until a 
weld pool the approximate width of the quasi-steady state pool formed and then input 
in a 4.5mm hemisphere. This heat input scheme caused cooling rates about 5 to 7 times 
higher than for the constant hemispherical input. This occurred because the heat input 
when concentrated on the surface caused temperatures in excess of 4000 K. When the 
heat input was shifted to the hemispherical distribution, the sudden diffusion of the heat 
input resulted in higher cooling rates. To alleviate this problem, the code was altered to 
input the heat on the surface until the weld pool width exceeded Imm, then input the 
heat in an ellipsoid of 4.5mm radius and 1.5mm depth. As the weld pool width at the 
surface increased, the depth of the heat input was also increased until the depth was 3.5 
mm. Fora heat input of 1950 watts, the increases in depth of the heat input distribution 
occurred at 0.08, 0.30, and 2.22 seconds after initiation of the arc. A comparative plot 
of the two cases above is provided in Figure 6. The parameters used were 1950 watts 
and 4mm per sec for both plots. The first plot is for a 15 second run with no cooldown; 
the second for seven seconds of welding followed by eight seconds with no heat input. 
In the first case, cooling rates of nearly 1400°C per second were predicted near the point 
of arc initiation. In the second case, the predicted cooling rates were less than 400°C 
per second at the same point. 

Comparing the results of the two cases above with those of Ule [Ref. 12: pp. 44-45], 
(where a surface input of heat and a directional thermal conductivity in the weld pool 
were utilized), it became apparent that the cooling rates determined by the model during 
the initiation of welding were heavily dependent on the heat distribution chosen. To 
adequately model this part of the weld cycle without taking convection in the weld pool 
as well as other important phenomena into consideration would require an excessive 


computational effort using a pure conduction model. The derivation of an empirical 


relationship that could be used with minimal error would require in-depth study of the 
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shapes of weld pools from start-up to quasi-steady state for the range of materials and 
process parameters expected to be encountered. Nearly the same argument can be ap- 
plied to the cooling rates at arc shutdown. 

The use of the codes to correctly determine the actual cooling rates during start-up 
and shutdown without some accounting for convection and other effects does not appear 


to be viable. The quasi-steady state cooling rates determined by the model do not suffer 


from the same limitations. 





HI. EXPERIMENTAL APPARATUS 


A. WELDING MACHINE 

The welding machine utilized for all welds performed during experimentation was a 
Miller DC Welding Power Source, Model SR600/SCMI1A with Electroslope 3. A trav- 
erse was included in the setup which allowed for three-dimensional positioning of the 
electrode. The forward and reverse directions were controlled by a motor. The trans- 
verse movement was manual and was only utilized during initial positioning of the torch. 
The vertical movement was also motor controlled, but not energized during welding. 
The welding machine was capable of providing direct current up to 675 amperes as well 
as a high frequency alternating current used to start the arc. The shielding gas used was 
pure Argon at flow rates of 15 to 20 cubic feet per hour. 


B. PLATE INSTRUMENTATION 

The seven- by twelve-inch, one inch thick HY-80 plates were instrumented with up 
to six platinum-platinum rhodium (type S) thermocouples, three mils in diameter. The 
thermocouple diameter was chosen as small as practicable to reduce response time. The 
placement of the thermocouples was as close to the arc centerline as the arc temperature 
would allow. The distance to the arc of the nearest thermocouple was 12 mm for the 
baseline runs and 10 mm for the surface flaw runs. The thermocouples were spaced at 
3 mm intervals. A diagram of the placement of the thermocouples is shown in 
Figure 7. The locations were marked by indenting the plates slightly with a punch. A 
high temperature, high thermal conductivity epoxy, CERAMABOND 569, was used to 
glue the thermocouples in place. The leads were thread through small ceramic capillaries 
which provided support for the fragile wires. A wooden support was manufactured for 
the capillary tubes and the connectors. This arrangement proved very durable and nu- 


merous runs were performed on the same plate. 


C. DATA ACQUISITION SYSTEM 

The thermocouples were connected to a high speed voltmeter with automatic com- 
pensation in an HP3852A Data Acquisition'Control Unit. The HP3852A was controlled 
by an HP9000 series 300 computer with a 9153C 20Mb hard disk and disk drive using 


Basic version 4.0. The program used to sample the data was existing and had been used 


extensively for welding laboratory experiments in support of a graduate course. 





Arc centerline . . 


> Arc centerline 





Figure 7. Thermocouple locations: (a) Baseline runs: offset 1l2ru, 3001 spacing 


(b) Surface flaw run: offset 1Or1, 3.0001 spacing 


The six thermocouple channels were scanned at 20 millisecond intervals for a period 
of about 40 seconds. The temperature versus time plots were observed on the terminal 
screen and also plotted for comparison with the computational model. A block diagram - 


of the welding and data acquisition svsterhs is provided in Figure $. 
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Figure 8. Welding and data acquisition systems 


D. ARC EFFICIENCY DETERMINATION 

The arc efficiency was measured by evaporating FC-72 when the just welded samole 
was placed in a tempering beaker of the fluid held near its saturation temperature. 
FC-72 is an inert dielectric liquid with a saturation temperature of about 56°C at one 
atmosphere pressure. It is often used in electronic cooling applications. The fluid was 


maintained at near saturation temperature by a heated recirculation bath with a flow 





rate of about 1.5 gallons per minute of distilled water. The bath was a NESLAB 
Endocal model RTE-5B. The energy to change the phase of the mass of fluid lost plus 
the energy to raise the sample to the saturation temperature of the fluid was equated 
with the heat input to the metal. The energy used by the machine was calculated from 
voltage, current and time measurements during the sample weld. The amount of fluid 
lost due to steady evaporation during the course of the experiment was measured with 
a balance with a least count of one gram. The efficiency was then calculated by dividing 
the average power input by the machine indication of power. A block diagram of the 


system used for the efficiency measurement is shown in Figure 9. 


E. FUSION ZONE COMPARISON 

The size and shape of the fusion zone during the transient from arc initiation to 
quasi-steady state was determined by welding on the sample, cutting the bead at speci- 
fied intervals, preparing the samples and photographing the prepared surfaces for 
measurements. The large plates were cut into sections small enough to fit into the 
sample cutter using a power hacksaw, then cut to the final size for analysis. The samples 
were prepared bv sanding on 240, 320, 400 and 600 grit emery boards followed bv pol- 
ishing with 0.05 alumina and etching for about one minute with a 2% Nitol solution. 


The fusion zones Were photographed using Polaroid film through an optical microscope. 





C. FUSION ZONE COMPARISON 

The model under-predicted the size of the fusion zones during start-up as shown in 
Figures 12 through 15. Results of measurements of the actual fusion zone at positions 
corresponding to various times after initiation of the arc are represented in comparison 
to the model predictions in Table 1. Additionally, the tabular data were plotted and 
appear in Figure 15. It was noted that initially the modeled fusion zone was wider than 
the actual fusion zone. This result would imply that in the model the heat may have 
been input over too wide an area at the beginning of the weld. As the weld progressed 
the actual and modeled results followed the same trend on the plot. One difference ob- 
served was the model fusion zone reached quasi-steady state at about three seconds 
while the actual fusion zone reached a size and shape about which it oscillated after 


about five seconds. 


Table 1. FUSION ZONE COMPARISON: Area ratio is the fusion zone area of 


the model divided by that of the sample. 
Area 
Ratio 


Sample Model 
Time Measurements Predictions 


(mm) (mm) (mm) (mm) 
L#s | sas |, zap, 0 es oe 
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DETERMINATION OF POWER INPUT INTO SAMPLE 
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Figure 10. Welding process efficiency 


plot. The results were therefore consistent with the assumptions of Rosenthal’s analytic 
solution and the structure of the model. 
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IV. EXPERIMENTAL RESULTS 
A. DETERMINATION OF POWER INPUT TO SAMPLE 


The efficiency determination was made by two different methods. The first method 
was boiling FC-72 and measuring the mass of fluid lost by evaporation when the just 
welded sample was placed in the fluid. A second method was used because the results 
of the first were not consistent. The second method was matching the temperature-time 
plots obtained from the thermocouple output during baseline runs to the Rosenthal 
solution for the same thermophysical properties and arc velocities. Since the measure- 
ments used for the comparison were well away from the fusion zone, the Rosenthal 
model was an appropriate first approximation. The results are shown in Figure 10. 

The inconsistent results of the phase change method may be explained in part by the 
relative magnitudes of the heat losses from the samples. The point furthest to the right 
in Figure 10 was obtained for about 14 seconds of welding, while the other two points 
were obtained from welds in excess of 35 seconds. For the latter two cases, the small 
sample size (2 inches by 5 inches by | inch thick) may have resulted in more conductive 
losses through the welding machine clamp as well as larger convective and radiative 
losses during welding and when the samples were moved to the tempering beaker. 

The Rosenthal equation matching was more consistent. An explanation for the 
deviations from a linear response is the uncertainties in the experimental method which 


are discussed at the end of this chapter. . 


B. TEMPERATURE-TIME PLOT COMPARISONS 

Several welds at various power levels were made and temperature-time plots ob- 
tained from the thermocouple output. The Rosenthal solution was utilized to match the 
actual power input to the metal with the temperature-time plots. When the START2 
program was run for the same parameters, the temperature-time plot did not match the 
experiment. In every case the model predicted much lower temperatures, though the 
shapes of the plots were similar. One source of the differences was that the Rosenthal 
solution does not include losses due to convection and radiation from the material as 
well as losses from the arc. Because the model applied the power directly to the sample 
and accounted for convective and radiative losses, a higher power than predicted by the 


Rosenthal solution would have been required to produce the same temperature-time 
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Figure 9. Efficiency determination system 
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Figure I]. Fusion zone comparison at 0.5 seconds: (a) Photograph at 16X mag- 
nification (b) Model prediction at 16X magnification, scale in mun 
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Figure 12. Fusion zone comparison at 1.5 seconds: (a) Photograph at 16X mag- 
nification (b) Model prediction at 16X magnification, scale in nun 
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Figure 13. Fusion zone comparison at 2.0 seconds: (a) Photograph at 16X mag- 
nification (b) Model prediction at 16X magnification, scale in mun 
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Figure 14. Fusion zone comparison at 3.0 seconds: (a) Photograph at 16X mag- 


nification (b) Model prediction at 16X magnification, scale in sun 
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Figure 15. 
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D. SURFACE FLAW COMPARISON 

The temperature-time plots obtained from the thermocouple output for the surface 
flaw experiment were not as graphic in their representation of the flaw as the model. 
The thermocouple output for the surface flaw experiment is plotted in Figure 16. The 
model] has the ability to show the temperature variations on the surface between the arc 
and the flaw. The limitations imposed by the use of surface mounted thermocouples 
precluded the same measurements in the experiment. The distortion of the thermal 
profile was evident, however, despite the noise in the thermocouple output. From this 
experiment it follows that a non-contact sensor capable of measuring the temperatures 
on the surface quite close to the arc, with a good sensitivity, would be more useful in 


validating the model capabilities in this area. 


E. EXPERIMENTAL UNCERTAINTIES 
Various components of the overall uncertainty are discussed in the following for the 
two phases of the experimental program. 
1. Efficiency determination 
The uncertainties in the procedure used to determine the power input to the 
sample during Welding due to measurements only occurred in the following parameters: 
time, current, voltage, and mass of the fluid. The uncertainties for each parameter were: 
time 0.5 seconds 
current 10 amps 
voltage 0.5 volts 


mass 0.5 grams 


The convective, radiative and conductive losses were unknown. 

For the Rosenthal equation matching there was approximately a one millimeter 
uncertainty in the location of the thermocouples relative to the arc due to the error in 
marking the plate for thermocouple locations and to deviations of the arc from its in- 
tended path caused by slop in the gears in the traverse. An additional source of error 
was in the thermocouple output itself. There was approximately 5 to 10 K variation in 
the output as shown in Figure 17. The plot was representative of the thermocouple 
output on all channels with the welding machine off and the sample at ambient tem- 


perature for several davs. The traverse speed was calibrated by timing the travel of the 


electrode over a distance of about 30 inches for each speed setting of the contro]. The 
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Figure 16. Surface flaw temperature-time plots: Flaw of 4num diameter and depth 
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response was linear but some error may have been encountered due to the coarse set- 
tings on the speed control. 
2. Fusion zone comparison 

The fusion zone comparison samples were cut with a blade of about 2mm width, 
which corresponded to about 0.5 seconds of elapsed time at a torch velocity of 4.56 
mm per second. Each sample was cut to about 4.5m thickness. The uncertainty in the 
measurements of the thickness was 0.5 mm. The uncertainty in the measurements to 
determine the exact location of the arc at a given time was 0.5mm. Combining these 
results, the actual time of the sample cut relative to the initiation of the arc had an un- 
certainty of about 0.25 seconds. The measurements of the fusion zone sizes on the 
photographs were highly accurate. The uncertainty of the measurements on the 16X 


magnified photographs was .03 inches. When correlated to the actual size of the sample, 


the uncertainty is fess than 0.05:77. 
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Figure 17. Noise in the thermocouple output 
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Vv. CONCLUSION AND RECOMMENDATIONS 
A. CONCLUSION 


The computational model written by Ule has been validated and shown to be useful 
in predicting the temperature distribution in a solid during autogenous welding. The 
ability to accurately predict the cooling rates during the weld cycle requires further de- 
velopment. The model has shown adaptability to different methods of heat input, which 
allow for study of different welding processes. The detectability of flaws in the material 
has been shown to be dependent on flaw location and composition. The need for an 


accurate non-contact surface temperature sensing system has been established. 


B. RECOMMENDATIONS 
The current study has accomplished several objectives, yet further study and re- 
search will be required to determine the full capabilities and usefulness of the computa- 
tional model. 
I. Cooling Rate Control 
Additional studies of the transient development of the weld pool during start-up 
and the solidification process during shutdown are required to determine the modifica- 
tions necessary in the heat input distribution for the model to adequately predict the 
cooling rates at those points in the process. Through comparison of observed fusion 
zone sizes and shapes with model predictions, such a modification could be developed 
in an iterative process. 
2. Heat Input Distribution 
The double elliptical heat distribution should be studied and modified in order 
to obtain the proper major and minor axes of the ellipses and ratio of energy input in 
front of and behind the arc according to the process and material in question. The re- 
sulting information would allow the code to be used for the entire range of autogenous 
processes. 
To minimize the complexity and subsequent run time of the code, convection in 
the weld pool should continue to be accounted for by use of a fictitiously high 


conductivity of the liquid phase and tke arc digging should continue to be accounted for 


by the proper choice of heat input distribution. 





3. Flaw Detectability 
The detectability of flaws for various materials and process variables should be 
determined to assist in the preliminary design of a non-contact surface temperature 
sensing system by providing an estimate of the required sensitivity. 
4. Temperature Sensing System 
The optimal non-contact surface temperature sensing system should be deter- 
mined. A laser vision system has shown a capability of viewing the weld pool of a sta- 
tionary arc. A qualitative determination of the shape and size of the fusion zone on the 
surface of the material for moving arcs should be made. This information could be ap- 


plied to the model to assist in determining the heat input distribution shape, especially 


during start-up. 





APPENDIX A. FLOW CHART 


A chart of the basic program flow was constructed to aid in the understanding of 


the code structure and execution. . 





Dimension zone arrays and summing arrays 


Declare variables 





Set welding parameters: voltage. amperage. efficiency, velocity and time 


Set material parameters: [Emissivity, density, specific heat capacity, hcat transfer 
cocflicient and Boltzmann’s constant 


Set arc heat input dimesions, fine zone grid size and center are on fine zone in x- . 
direction 


Determine thermal conductivity and thermal diffusivity for material 
t grid sizes of all grids 


Open output files 


Set imitial values of counters and pointers 















Initialize fine zone enthalpy array 





Calculate input power in watts and input energy in joules for each time step 











Begin time step loop 


Increment time 
Calculate distance traveled 

Step time loop counter 

Position arc 

Calculate weighting factor 

‘\dd heat to fine zone in specified manner 
Convert enthalpy to temperature in fine zone 

Call finite difference subroutines; calculate T + K 














1/2 
Call finite diflerence subroutines; calculate T + K2/2 
Call finite difference subroutines; calculate T + K3 
Call finite difference subroutines: calculate Runge-Kutta sum 
Check if fine grid stepping required; if so, step fine grid 
Check if medium grid stepping required; ifso, step medium grid 
Output solution files at specified time intervals 


End time step loop 


Output final results 






















Allow for option of restarting the code 


User delined functions: Temperature as function of enthalpy, enthalpy as a function 
of temperature and thermal conductivity as a function of temperature 





Figure [8 Program Flow Chart 





APPENDIX B. MODIFIED PROGRAMS ; . 


The portions of codes that were changed are provided along with the corresponding 
portions of Ule’s codes for comparison. A list of the additional variables used in each 


new section of code is also provided. 


A. BASIC MODIFICATIONS 
1. START2 Program Modified Sections 


SIGMA = 5.67E-8 
EPSIL = 0. 82 
RHO = 7890 
CPSP = 536 
HF = 25.0 
HM = 10.0 
HC = 10.0 
C 
Al = 4.0 
Bl = 4.0 
C1 = 4.0 
SPACE = 1.0 
XARC = 14.0 
C 
KAPPA = FK(300. 0) : 
ALPHA = KAPPA/(RHO*CPSP) 
C 
DELXF = SPACE . 
DELXM = 3.0 * SPACE 
DELXC = 3.0 * DELXM 
DELYF = SPACE 
DELYM = 3.0 * SPACE 
DELYC = 3,0 * DELYM 
DELZF = SPACE 
DELZM = 3.0 * SPACE 
DELZC = 3.0 * DELZM 


C BOUNDARY CONDITION COEFFICIENTS 


FO( 1)=ALPHA*DELT/(DELXC*DELYC*1. 0E-3**2) 
B1(1)=(2. 0/3. 0)*HC*DELXC*1. OE-3/KAPPA 


FOC 2)=ALPHA*DELT/( DELXM*DELYM*1. OE-3**2) 
BI(2)=2. O*HM*DELXM/KAPPA 


FO(3)=DELT/(DELXF*1. 0OE-3)**2 
B1I(3)=EPSIL*SIGMA*2, 0*DELT/DELZF 
BI(4)=HF*2. 0*DELT/DELZF 


C CALCULATE THE RUNGE-KUTTA APPROXIMATION 





DO 10 M=1,NDIV 
TIME=FLOAT(NINT( 100. *(TIME+DELT) ))/100. 
DIS=TIME*VEL 
. N=N+1 


C POSITION HEAT SOURCE, CALCULATE VOLUME WEIGHTING FACTOR 'ENERGY' 


YARC=DIS+73-NB*9 -NC*3 
ENERGY = 0.0 
I1 = INT (A1/SPACE) 
XLIML = XARC - I1-1 
XLIMU = XARC + 1141 
J1 = INT(B1/SPACE) 
YLIML = INT(YARC) - J1 - 2 
YLIMU = INT (YARC) + J1 + 2 
K1 = INT(C1/SPACE) + 1 
DO 300 I = XLIML,XLIMU 
XE = ((1-XARC)/A1)**2 
DO 200 J = YLIML,YLIMU 
YE = ((J-YARC)/B1)**2 
DO 100 K = 1,Kl 
ZE = ((K-1)/C1)**2 
EXPON = XE+YE+ZE 
IF(EXPON. GT. 1.0) GO TO 200 
IF(K.EQ.1) THEN 
VOL = 0.5 * (SPACE*1. 0OE-3)**3 


» ELSE 
VOL = (SPACE*1. 0E-3)**3 
ENDIF 
ENERGY = ENERGY + VOL * EXP(-3.0 * EXPON) 
100 CONTINUE 


200 CONTINUE 
_300 CONTINUE 


C ADD THE HEAT FROM THE ARC HERE USING A SIMILAR DO LOOP 
C CONSTRUCTION AS ABOVE! 


DO 600 I = XLIML,XLIMU 
XE = ((1-XARC)/A1)**2 
DO 500 J = YLIML,YLINU 
YE = ((J-YARC) /B1)**2 
DO 400 K = 1,K1 
ZE = ((K-1)/C1)**2 
EXPON = XE+YE+ZE 
IF(EXPON. GT. 1.0) GO TO 500 
C(1I,J,K)=C(1,J,K)+EXP(-3. O*EXPON)*Q/ENERGY 
400 CONTINUE 
- 500 CONTINUE 
600 CONTINUE 


4] 











2. Variable List 
The following variables were added to all programs to allow easier input of dif- 
ferent heat input schemes and to allow for the study of different materials. The grid size . . 


parameters were also coded for future use if it was desired to change grid size. All vari- 


ables are given in standard SI units unless otherwise noted. : 
SIGMA Boltzmann's constant 
EPSIL Emmissivity of the surface 
RHO Density of the material 
CPSP Specific heat capacity 
HF Heat transfer coefficient for the fine zone surfaces 
HM Heat transfer coefficient for the medium zone surfaces 
HC Heat transfer coefficient for the coarse zone surfaces 
Al Radius of heat input in x-direction in millimeters 
Bi Radius of heat input in positive y-direction in millimeters 
B2 Radius of heat input in negative v-direction in millimeters 
Cl Radius of heat input in z-direction in millimeters 
SPACE Node spacing in the fine zone in millimeters 
XARC Position of the arc in the fine zcne in x-direction ° 
KAPPA Thermal conductivity 
ALPHA Thermal diffusivity 
DELXF Control] volume length in x-direction in fine zone 
DELXM Control] volume length in x-direction in medium zone 
DELXC Control] volume length in x-direction in coarse zone 
DELYF Contro] volume length in y-direction in fine zone 
DELYM Control] volume length in y-direction in medium zone 
DELYC Control] volume length in y-direction in coarse zone 
DELZF Control volume length in z-direction in fine zone 
DELZM Control volume length in z-direction in medium zone 
DELZC Control volume length in z-direction in coarse zone 
FO(1) Fourier number in the coarse zone 
BI()) Biot number in the coarse zone - 
FO(2) Fourier number in the medium zone 


Bl(2) Biot number in the medium zone 





FO(3) Fourier number in the fine zone 


BI(3) Biot number in the fine zone (for radiation) 

BI(4) Biot number in the fine zone (for convection) 

ENERGY Summing variable for the total energy input in a time step 

I] Converts radius of arc input to nodal value for x-direction 
XLIML Lower limit of x-direction input of heat 

XLIMU Upper limit of x-direction input of heat 

J1 Converts radius of arc input to nodal value for y-direction 

YLIML Lower limit of y-direction input of heat 

YLIMU Upper limit of y-direction input of heat 

K1 Converts radius of arc input to nodal value for z-direction 

XE Contribution to exponent in Gaussian distribution from x-direction 
YE Contribution to exponent in Gaussian distribution from y-direction 
ZE Contribution to exponent in Gaussian distribution from z-direction 
EXPON Value of exponent in the Gaussian distribution expression 

VOL Volume of the node being evaluated in cubic meters 


3. Original WELD Program Sections 
The following are sections of program from the original version written and used 
by Ule [Ref 12: pp. 85-86]. 


DATA XZ/. 1433987,. 2870068, 
%,4352117,.5,.4352117,. 2870068,. 1433987,. 2184978,. 4373148, 
*, 6631358,. 76118543,. 6631358,.4373148,. 2184978,. 096619,. 1933793, 
*, 2932366,. 3368896,. 2932366,. 1933793,.096619/ 

DATA A,B/8046*300. 0/,C/5832*1. 14237E8/,ASUM, BSUM/8046*0. 0/ 
* , AOUT , BOUT /8046*0. 0/ ,CSUM/5832*0. / ,COUT/5832*0. / 


C INITIAL DATA BLOCK FOR STARTING A PROBLEM, SET PROBLEM LENGTH AND 
C ARC PARAMETERS 


FINI = 10.0 
VOLT=30. 
AMP=265. 
EFF=. 32 
VEL=4. 


C OPEN THE OUTPUT FILES 
STATUS='NEW' , FORM='UNFORMATTED' ) 
,STATUS='NEW , FORM='UNFORMATTED' ) 


OPEN(1,FILE='SURF' 


OPEN(2,FILE='FINAL' 
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OPEN(3,FILE='CUT' 


STATUS='NEW' , FORM='UNFORMATTED' ) 
OPEN(4,FILE='HIST 


 STATUS=' NEW , FORM='UNFORMATTED' ) 
C THE INITIAL CONDITIONS: NUMBER OF DIVISIONS, TIME STEP, TIME . . | 
C AND ARC LOCATION, ETC. 


NDIV=FINI*100 
DELT=. 01 . 
TIME=0. 

STEP=3. 

OUT=. 49 

N=0 

BSTEP=0 

CSTEP=0 

NB=3 

NC=10 

QDENSE=7. 134346E-9 
TINF=300. 0 


C WELD PARAMETERS 


QDOT=EFF*VOLT*AMP 
Q=QDOT*DELT 


C BOUNDARY CONDITION COEFFICIENTS 


FO(1)=DELT*. 1636 

FO(2)=DELT*1. 4722 

FO(3)=DELT*1000000. ’ 
BI(1)=. 001132 

BI(2)=. 001132 

BI(3)=DELT*. 00009299 

BI(4)=DELT*50000. ’ 


C CALCULATE THE RUNGE-KUTTA APPROXIMATION 
DO 10 M=1,NDIV 
TIME=TIME+DELT 
DIS=TIME*VEL 
N=N+1 
C POSITION HEAT SOURCE AND CALCULATE VOLUME WEIGHTING FACTOR SUM! 


YARC=VEL*TIME+73-NB*9-NC*3 


SUM=0. 
DO 1 J=7,23 
c IF ((J-YARC).GT.(0.0)) THEN LINES ALLOWED SHAPING THE ARC 
SUM=SUM+QDENSE/4. *EXP( -. 10625%((J-YARC)**2) ) 
C ELSE 
C SUM=SUM+QDENSE/10. *EXP( -. 017*((J-YARC)**2)) 
C ENDIF 
1 CONTINUE ; 
SUM=SUM/Q 


C ADD THE HEAT FROM THE ARC 





DO 2 J=7,23 
IF ((J-YARC).GT.(0.0)) THEN + LINES ALLOWED SHAPING THE ARC 
Y=. 25*EXP(-. 10625*((J-YARC)**2)) 
ELSE 
Y=. 1*EXP( -. 017*(( J-YARC)**2)) 
ENDIF 
DO 2 I=11,17 
: DO 2 K=1,3 
C(I,J,K)=C(1,J3,K)+Y*XZ(I-10,K)/SUM 
2 CONTINUE 
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